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Pressure-energy correlations in liquids. I. Results from computer simulations 

Nicholas P. Bailey/ '[jUlfR. Pedersen/ Nicoletta Gnan/ Thomas B. Sclir0der/ and Jeppe C. Dyre^ 

^DNRF Center "Glass and Time'', IMFUFA, Dept. of Sciences, 

Roskilde University, P.O. Box 260, DK-4000 Roskilde, Denmark 

(Dated: October 22, 2008) 

We show that a number of model liquids at fixed volume exhibit strong correlations between equilibrium 
fluctuations of the configurational parts of (instantaneous) pressure and energy. We present detailed results 
for thirteen systems, showing in which systems these correlations are significant. These include Lennard- Jones 
liquids (both single- and two-component) and several other simple liquids, but not hydrogen-bonding liquids like 
methanol and water, nor the Dzugutov liquid which has significant contributions to pressure at the second nearest 
neighbor distance. The pressure-energy correlations, which for the Lennard- Jones case are shown to also be 
present in the crystal and glass phases, reflect an effective inverse power-law potential dominating fluctuations, 
even at zero and slightly negative pressure. An exception to the inverse-power law explanation is a liquid with 
hard- sphere repulsion and a square- well attractive part, where a strong correlation is observed, but only after 
time- averaging. The companion paper [arXiv: 0807. 0551] gives a thorough analysis of the correlations, with a 
focus on the Lennard- Jones liquid, and a discussion of some experimental and theoretical consequences. 



I. INTRODUCTION 

Physicists are familiar with the idea of thermal fluctuations 
in equilibrium. They also know how to extract useful infor- 
mation from them, using linear response theory.^ ^^^ These 
methods started with Einstein's observation that the specific 
heat in the canonical ensemble is determined by the magni- 
tude of energy fluctuations. In any thermodynamic system 
some variables are fixed, and some fluctuate. The magnitude 
of the variances of the latter, and of the their mutual covari- 
ance, determine the thermodynamic "response" parameters.^ 
For example, in the canonical (NVT) ensemble, pressure p and 
energy E fluctuate; the magnitude of pressure fluctuations is 

related to the isothermal bulk modulus Kt = —V {^\ , 

that of the energy fluctuations to the specific heat at constant 
volume cy = ^ (ff )y' while the covariance (ApAE) is 

related"^ to the thermal pressure coefficient Py = (^j .If 

the latter is non-zero, it implies a degree of correlation be- 
tween pressure and energy fluctuations. There is no obvious 
reason to suspect any particularly strong correlation, and to 
the best of our knowledge none has ever been reported. But 
in the course of investigating the physics of highly viscous 
liquids by computer simulation, we noted strong correlations 
between pressure and energy equilibrium fluctuations in sev- 
eral model liquids, also in the high temperature, low- viscosity 
state. These included the most studied of all computer liquids, 
the Lennard- Jones system. Surprisingly, these strong correla- 
tions survive crystallization, and they are also present in the 
glass phase. "Strong" here and henceforth means a correla- 
tion coefficient of order 0.9 or larger. In this paper we ex- 
amine several model liquids and detail which systems exhibit 
strong correlations and which do not. In the companion papei - 
(referred to as Paper II) we present a detailed analysis of the 
correlations for the single-component Lennard- Jones system, 
and discuss some consequences. 

Specifically, the fluctuations which are in many cases 
strongly correlated are those of the configurational parts of 
pressure and energy (the parts in addition to the ideal gas 



terms). The (instantaneous) pressure p and energy E have 
contributions both from particle momenta and positions: 



p = NkBT{p,, . . . , pn)/V + W{t,, . . . , tn)/V 

E = K(pi, . . . , Piv) + U{ri, . . . , r^), (1) 

where K and U and the kinetic and potential energies, re- 
spectively. Here T(pi, . . . , p^v) is the "kinetic temperature"!^ 
proportional to the kinetic energy per particle. The config- 
urational contribution to pressure is the virial W, which is 
defined^ by 



W 



\^r,-Vr^U 



(2) 



where r^ is the position of the ith particle. Note that W has 
dimension energy. For a pair interaction we have 



Upair = ^V{rij) 



(3) 



Kj 



where r^j is the distance between particles i and j and v{r) 
is the pair potential. The expression for the virial (Eq. (|2])) 
become^ 



V^pair = -i^^rijv'iuj) = --^w{rij) (4) 



i<j 



where for convenience we define 



i<j 



w{r) = rv' {r). 



(5) 



Fig. [T] (a) shows normalized instantaneous values of p and 
E, shifted and scaled to have zero mean and unit variance, as a 
function of time for the standard single-component Lennard- 
Jones (SCLJ) liquid, while Fig.[T](b) shows the corresponding 
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FIG. 1: (Color online) Equilibrium fluctuations of (a) pressure p 
and energy E and (b) virial W and potential energy [/, in a single- 
component Lennard- Jones system simulated in the NVT ensemble at 
p = 34.6 mol/1 and T = 80K (Argon units). The time-averaged 
pressure was close to zero (1.5 MPa). The correlation coefficient R 
between W and U is 0.94, whereas the correlation coefficient is only 
0.70 between p and E. Correlation coefficients were calculated over 
the total simulation time (10 ns). 



FIG. 2: (Color online) (a) Scatter-plot of instantaneous virial W and 
potential energy U from the simulation of Fig.[T] The dashed line is 
a guide to the eye, with a slope determined by the ratio of standard 
deviations of W and U (Eq. (|7|). (b) Example of a system with 
almost no correlation between W and U: TIP5P water at T=12.5°C, 
density 1007.58 kg/m^ (NVT). This system has Coulomb, in addition 
to Lennard- Jones, interactions. 



fluctuations of W and U. We quantify the degree of correla- 
tion by the standard correlation coefficient R, defined by 



R 



(AWAU) 



./mmvmm' 



(6) 



Here angle brackets () denote thermal averages while A de- 
notes deviation from the average value of the given quan- 
tity. The correlation coefficient is ensemble-dependent, but 
our main focus — the R ^ 1 limit-is not. Most of the sim- 
ulations reported below were carried out in the NVT ensem- 
ble. Another important characteristic quantity is the "slope" 
7, which we define as the ratio of standard deviations: 






(7) 



Considering the "total" quantities, p and E, (Fig. [TJa)) there 
is some correlation; the correlation coefficient 0.70. For the 
configurational parts, W and [/, on the other hand (Fig.[TJb)), 
the degree of correlation is much higher, R = 0.94 in this 
case. Another way to exhibit the correlation is a scatter-plot 
of W against U, as shown in Fig.[2ja). 

Is this correlation surprising? Actually, there are some in- 
teratomic potentials for which there is a 100% correlation 
between virial and potential energy. If we have a pair po- 
tential of the form v{r) ex r~^, an inverse power-law, then 
w{r) = —nv{r) and Wpair = (^/3)t/pair holds exactly. In this 
case the correlation is 100% and 7 = n/3. 



^ 1 

OJ 1 

W 






1 n r 




Lennard- Jones 

4/3 (r/a)'^^ - 4/3 



Repulsive part 

Attractive part 



1 r^ 1.2 1.4 

Distance/o 



1.6 



FIG. 3: (Color online) Illustration of the "effective inverse power- 
law" chosen in this case to match the Lennard- Jones potential and its 
first two derivatives at the point r = a. The vertical line marks the 
division into the repulsive and attractive parts of the Lennard- Jones 
potential. 



Conversely, suppose a system is known to be governed by 
a pair potential, and that there is 100% correlation between 
W and U. We can write both U and W at any given time t 
as integrals over the instantaneous radial distribution function 
defined^ as 






from which 



U{t) 



N 



p / dr 47Tr'^g{r^t)v{r) 



(9) 



and 



W(t) = p dr47rr^g(r,t)w(r). (10) 

6 Jo 

Here the factor of ^ is to avoid double-counting, and p = 
N/V is the number density. 100% correlation means that 
W{t) = 7t/(t) holds for arbitrary ^(r, t) (a possible addi- 
tive constant could be absorbed into the definition of U). In 
particular we could consider g{r^t) = S{r — ro)P Substi- 
tuting this into the above expressions, the integrals go away 
and we find w{ro) = —3'yv{ro). Since r^ was arbitrary. 



v'{r) 



-3'yv{r)/r, which has the solution v{r) ex 



-37^ 

This connection between an inverse power-law potential and 
perfect correlations suggests that strong correlations can be 



attributed to an effective inverse power-law potential, with ex- 
ponent given by three times the observed value of 7. This will 
be detailed in Paper II, which shows that while this expla- 
nation is basically correct, matters are somewhat more com- 
plicated than this. For instance, the fixed volume condition, 
under which the strong correlations are actually observed, im- 
poses certain constraints on ^(r, t). 

The celebrated Lennard- Jones potential is given^by 



VLj{r) = 4e 



(7)" -(7)' 



(11) 



One might think that in the case of the Lennard- Jones poten- 
tial the fluctuations are dominated by the repulsive r~^^ term, 
but this naive guess leads to a slope of four, rather than the 6.3 
seen in Fig.|2](a). Nevertheless the observed correlation, and 
the above mentioned association with inverse power-law po- 
tentials, suggest that an effective inverse power-law descrip- 
tion (involving short distances), with a more careful identifi- 
cation of the exponent, may apply. In fact, the presence of 
the second, attractive, term, increases the effective steepness 
of the repulsive part, thus increasing the slope of the correla- 
tion, or equivalently the effective inverse power-law exponent 
(Fig.|3]). Note the distinction between repulsive term and re- 
pulsive part of the potential: the latter is the region where 
v{r) has negative slope, thus the region r < rm (Vm being 
the distance where the pair potential has its minimum, 2^/^(7 
for vlj)- This region involves both the repulsive and attrac- 
tive terms (see Fig. [3] which also illustrates the approximation 
of the repulsive part by a power law with exponent 18). The 
same division was made by Weeks, Chandler, and Andersen 
(WCA) in their noted paper of 1971,^ in which they showed 
that the the thermodynamic and structural properties of the 
Lennard- Jones fluid were dominated by the repulsive part at 
high temperatures for all densities, and also at low tempera- 
tures for high densities. Ben-Amotz and Stell^ noted that the 
repulsive core of the Lennard- Jones potential may be approx- 
imated by an inverse power-law with n ~ 18-20. The approx- 
imation by an inverse power-law may be directly checked by 
computing the potential and virial with an inverse power-law 
potential for configurations drawn from actual simulations us- 
ing the Lennard- Jones potential. The agreement (apart from 
additive constants) is good, see Paper II. 

Consider now a system with different types of pair inter- 
actions, for example a binary Lennard- Jones system with A A, 
BB, and AB interactions, or a hydrogen-bonding system mod- 
elled via both Lennard- Jones and Coulomb interactions. We 
can write arbitrary deviations of U and W from their mean 
values, denoted AU and AW, as a sum over types (indexed 
by t; sums over pairs of a given type are implicitly under- 
stood): 



AU = ^AUu AW = ^AWt. 



(12) 



Now, supposing there is near-perfect correlation for the in- 
dividual terms with corresponding slopes 7^, we can rewrite 
Al^as 



/^W = Y,lt^Ut. 



(13) 



If the 7t are all more or less equal to a single value 7, then 
this can be factored out and we goi /SW :^ 7 At/. Thus 
the existence of different Lennard- Jones interactions in the 
same system does not destroy the correlation, since they have 
7t ^ 6. On the other hand the slope for Coulomb interaction, 
which as an inverse power-law has perfect W^ U correlations, 
is 1/3, so we cannot expect overall strong correlation in this 
case (Fig.[2](b)). Indeed such reasoning also accounts for the 
reduction of correlation when the total pressure and energy 
are considered: /S.E = A/7 + Ai^, while (for a large atomic 
system) VAp = 7A/7 + (2/3)Air. The fact that 7 is (for 
the Lennard- Jones potential) quite different from 2/3 implies 
that the p, E correlation is significantly weaker that of W^ U 
(Fig.[T]). Even in cases of unequal slopes, however, there can 
be circumstances under which one kind of term, and there- 
fore one slope, dominates the fluctuations. In this case strong 
correlations will be observed. Examples include the high- 
temperature limits of hydrogen-bonded liquids (section IIID[ ) 
and the time-averaged (total) energy and pressure in viscous 
liquids (Paper II). 

Some of the results detailed below were published previ- 
ously in Letter form;^ the aim of the present contribution is to 
make a comprehensive report covering more systems, while 
Paper II contains a detailed analysis and discusses applica- 
tions. In the following section, we describe the systems simu- 
lated. In section |lll| we present the results for all the systems 
investigated, in particular the degree of correlation (correla- 



tion coefficient R) and the slope. Section IV gives a summary. 



II. SIMULATED SYSTEMS 

A range of simulation methods, thermodynamic ensembles 
and computational codes were used. One reason for this was 
to eliminate the possibility that the strong correlations are 
an artifact of using a particular ensemble or code. In ad- 
dition, no one code can simulate the full range of systems 
presented. Most of the data we present are from molecular 
dynamics (MD) simulations, although some are from Monte 
Carlo^^ (MC) and event-driven^^ (ED) simulations. Most of 
the MD simulations (and of course all MC simulations), had 
fixed temperature (NVT), while some had fixed total en ergy 
(NVE). Three MD codes were used: Gromacs {GRO)P^ 
Asap (ASAP),^"^ and DigitalMaterial (DM).^^ Home-made 
(HM) codes were used for the MC and ED simulations. 

We now list the thirteen systems studied, giving each 
a code-name for future reference. The systems include 
monatomic systems interacting with pair potentials, binary 
atomic systems interacting with pair potentials, molecular sys- 
tems consisting of Lennard- Jones particles joined rigidly to- 
gether in a fixed configuration (here the Lennard- Jones inter- 
action models the van der Waals forces), molecular systems 
which have Coulomb as well as Lennard- Jones interactions, 
metallic systems with a many-body potential, and a binary 



system interacting with a discontinuous "square- well" poten- 
tial. Included with each system is a list specifying which 
simulation method(s), which ensemble(s), and which code(s) 
were used (semi-colons separate the method(s) from the en- 
semble(s) and the ensemble(s) from the code(s)). Details of 
the potentials are given in Appendix [A| 

CU Pure liquid Cu simulated using the many-body poten- 
tial derived from effective medium theory (EMT)P^Ell 
(MD; NVE; ASAP) 

DB Asynmietric "dumb-bell" molecules, ^^ consisting of two 
unlike Lennard- Jones spheres connected by a rigid 
bond; (MD; NVT; GRO) 

DZ The potential introduced by Dzuguto\^^ as a candidate 
for a monatomic glass-forming system. Its distinguish- 
ing feature is a peak mv{r) around 1.5 a, after which it 
decays exponentially to zero at a finite value of r; (MD; 

NVT, NVE; DM) 

EXP A system interacting with a pair potential with expo- 
nential repulsion and a van der Waals attraction; (MC; 
NVT; HM) 

KABLJ The Kob-Andersen binary Lennard- Jones liquidp^ 
(MD; NVT, NVE; GRO, DM) 

METH The Gromos^^ 3-site model for methanol; (MD; 
NVT; GRO) 

MGCU A model of the metallic alloy MggsCuis using an 
EMT-based potential;^^ (MD; NVE; ASAP) 

OTP A three-site model of the fragile glass-former Ortho- 
terphenyl (OTP);23(MD; NVT; GRO) 

SCLJ The standard single-component Lennard- Jones system 
with the interaction given in Eq. ([TT]); (MD, MC; NVT, 
NVE; GRO, DM) 

SPC/E The SPC/E model of water;^^ ^y^^. ^yj. qrq) 

SQW A binary model with a pair interaction consisting 
of an infinitely hard core and an attractive square 

TIP5P A five-site model for liquid water which reproduces 
the density anomaly.^^ (MD; NVT; GRO) 

TOL A 7-site united-atom model of toluene; (MD; NVT; 
GRO) 

The number of particles (atoms or molecules) was in the range 
500-2000. Particular simulation parameters (A^, p, T, dura- 
tion of simulation) are given when appropriate in the results 
section. 
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FIG. 4: (Color online) Scatter plots of the configurational parts of pressure and energy - virial versus potential energy - for several state points 
of the single-component Lennard- Jones liquid (NVT). Each oval represents simulations at one particular temperature and density where each 
data point marks instantaneous values of virial and potential energy. The dashed lines mark constant density paths with the highest density to 
the upper left (densities: 39.8 mol/1, 37.4 mol/1, 36.0 mol/1, 34.6 mol/1, 32.6 mol/1). State points on the dotted line have zero average pressure. 
The plot includes three crystallized samples (lower left corner), discussed at the end of section[lIIA|and, in more detail, in Paper II [reproduced 
fromRef.m. 



III. RESULTS 

A. The standard single-component Lennard- Jones system 

SCLJ is the system we have most completely investigated. 
VF, [/-plots are shown for a range of thermodynamic state 
points in Fig. |4] Here the ensemble was NVT with N=864, 
and each simulation consisted of a 10 ns run taken after 10ns 
of equilibration; for all SCLJ results so-called "Argon" units 
are used (a = 0.34 nm, e = 0.997 kJ/mol). Each elongated 
oval in Fig. |4] is a collection of W^ U pairs for a given state 
point. Varying temperature at fixed density moves the oval 
parallel to itself, following an almost straight line as indicated 
by the dashed lines. Different densities correspond to differ- 
ent lines, with almost the same slope. In a system with a pure 
inverse power-law interaction, the correlation would be exact, 
and moreover the data for all densities would fall on the same 
straight line (see the discussion immediately after Eq. ([5])). 
Our data, on the other hand, show a distinct dependence on 
volume, but for a given volume, because of the strong corre- 



lation, the variation in W is almost completely determined by 
that of U. 

Values of correlation coefficient R for the state points of Fig.|4] 
are listed in Table[l| along with the slope 7. In Fig. [5] we show 
the temperature dependence of both R and 7 for different den- 
sities. Lines have been drawn to indicate isochores and one 
isobar (p = 0). Note that when we talk of an isobar here, we 
mean a set of NVT ensembles with F, T chosen so that the 
thermal average of p takes on a given value, rather than fixed- 
pressure ensembles. This figure makes it clear that for fixed 
density, R increases as T increases, while it also increases 
with density for fixed temperature; the slope slowly decreases 
in these circumstances. In fact it eventually reaches four, the 
value expected for a pure r~^^ interaction (e.g., at p = 34.6 
mol/1, T = 1000 K, 7 = 4.61, see Ref. 9). This is consis- 
tent with the idea that the repulsive part, characterized by an 
effective inverse power-law, dominates the fluctuations: in- 
creasing either temperature or density increases the frequency 
of short-distance encounters while reducing the typical dis- 
tances of such encounters. On the other hand, along an isobar. 



TABLE I: Correlation coefficients R and effective slopes 7 for the 
single-component Lennard- Jones system (SCLJ) for the state points 
in Fig. |4] p is the thermally averaged pressure. The last five states 
were chosen to approximately follow the isobar p = 0. 
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these two effects work against each other, since as T increases, 
the density decreases. The density effect "wins" in this case, 
which is equivalent to a statement about the temperature and 
volume derivative of R: Our simulations imply that 
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where a^ = {dV/dT)p/V is the thermal expansivity at con- 
stant pressure and p is the particle density. This can be 
recast in terms of logarithmic derivatives (valid whenever 
{dR/dp)T > 0) as follows 
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FIG. 5: (Color online) Upper plot, correlation coefficient R for 
the SCLJ system as a function of temperature for several densities 
(NVT). This figure makes clear the different effects of density and 
temperature on R. Lower plot, effective slope 7 as a function of T. 
Simulations at temperatures higher than those shown here indicate 
that the slope slowly approaches the value four as T increases. This is 
to be expected because as collisions become harder, involving shorter 
distances, the effective inverse power-law exponent approaches the 
12 from the repulsive term of the Lennard- Jones potential. 
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(16) 



Thus what we observe in the simulations, namely that the 
correlation becomes stronger as temperature is reduced at 
fixed pressure, is a priori more to be expected when the ther- 
mal expansivity is large (since then the right hand side of 



Eq. ([16]) is large). This has particular relevance in the con- 
text of supercooled liquids, which we discuss in paper II, be- 
cause these are usually studied by lowering temperature at 
fixed pressure. On the other hand if the expansivity becomes 
small, as for example, when a liquid passes through the glass 
transition, the inequality ([16]) is a priori less likely to be sat- 
isfied. We have in fact observed this in a simulation of OTP: 
upon cooling through the (computer) glass transition, the cor- 
relation became weaker with further lowering of temperature 
at constant pressure. 

Remarkably, the correlation persists when the system has 
crystallized, as seen in the data for the highest density — the 
occurrence of the first-order phase transition can be inferred 
from the gap between the data for 90K and 11 OK, but the 
data fall on the same line above and below the transition. 
One would not expect the dynamical fluctuations of a crystal. 
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FIG. 6: (Color online) Scatter-plot of the W, U correlations for a 
perfect face-centered cubic (FCC) crystal of Lennard- Jones atoms 
at temperatures IK, 2K, 3K, 5K, lOK, 20K, 30K, 40K, 50K, 60K, 
70K and 80K, as well as for defective crystals (i.e., crystallized from 
the liquid) at temperatures 50K, 70K and 90K (NVT). The dashed 
line gives the best fit to the (barely visible) lowest-temperature data 
(T = IK). The inset shows the temperature dependence of R at very 
low temperatures. The crystalline case is examined in detail in Paper 
II, where we find that R does not converge to unity at T = 0, but 
rather to a value very close to unity. 



which are usually assumed to be well-described by a harmonic 
approximation, to resemble those of the high-temperature liq- 
uid. In fact for a one-dimensional crystal of particles interact- 
ing with a harmonic potential v{r) = ^/c(r — r^)^ it is easy to 
show (Paper II) that there is a negative correlation with slope 
equal to -2/3. To investigate whether the harmonic approx- 
imation ever becomes relevant for the correlations, we pre- 
pared a perfect FCC crystal of SCLJ particles at zero temper- 
ature and simulated it at increasing temperatures, from 0.02K 
to 90K in Argon units, along a constant density path. The re- 
sults are shown in Fig.|6] Clearly the correlation is maintained 
right down to zero temperature. The harmonic approximation 
is therefore useless for dealing with the pressure fluctuations 
even as T ^ 0, because the slope is far from -2/3. The reason 
for this is that the dominant contribution to the virial fluctua- 
tions comes from the third order term, as shown in Paper II. 



B. A case with little correlation: the Dzugutov system 



FIG. 7: (Color online) A plot of the Dzugutov pair potential, with 
the Lennard- Jones potential (shifted by a constant) shown for com- 
parison. 
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FIG. 8: (Color online) Scatter-plot of W, U correlations for the 
Dzugutov system at density 0.88 and temperatures 0.65 and 0.70 
(NVE). The dashed lines indicate the best-fit line using linear regres- 
sion. These are consistent with the temperature dependence of the 
mean values of (W) and ([/), as they should be (see appendix [b]), 
but they clearly do not represent the direction of greatest variance. 
The full lines have slopes equal to the ratio of standard deviations of 
the two quantities (Eq. ([T])). The correlation coefficient is 0.585 and 
0.604 for T = 0.65 and T = 0.70, respectively. 



Before presenting data for all the systems studied, it is use- 
ful to see what it means for the correlation not to hold. In this 
subsection we consider the Dzugutov system,E^ whose poten- 
tial contains a peak at the second-neighbor distance, (Fig. [7] 
see appendix|A|for details) whose presence might be expected 
to interfere with the effectiveness of an inverse power-law 
description. In the next subsection we show how in a spe- 
cific model of water the lack of correlation can be explicitly 
seen to be the result of competing interactions. Fig. [8] shows 
W^ U plots for the Dzugutov system for two nearby tempera- 



tures at the same density. The ovals are much less elongated 
than was the case for SCLJ, indicating a significantly weaker 
correlation — the correlation coefficients here are 0.585 and 
0.604, respectively. In paper II it is shown explicitly that the 
weak correlation is due to contributions arising from the sec- 
ond peak. Note that the major axes of the ovals are not aligned 
with the line joining the state points, given by the mean val- 
ues of W and U, here identifiable as the intersection of the 
dashed and straight lines. On the other hand, the lines of best 



fit from linear regression, indicated by the dashed fines in each 
case, do coincide with the fine connecting state points. This 
holds generally, a fact which follows from statistical mechan- 
ics (appendix [B]). The interesting thing is rather that the major 
axes point in different directions, whereas in the SCLJ case 
they are also aligned with the state-point line. The linear- 
regression slope, being equal to (A/7Aiy)/((A/7)^), treats 
W and U in an asymmetric manner by involving ((A/7)^), 
but not ((AVK)^). This is because a particular choice of in- 
dependent and dependent variables is made. If instead we 
plotted U against W, we would expect the slope to be sim- 
ply the inverse of the slope in the W^ U plot, but in fact the 
new slope is {/SU /SW) / {{/SW)'^) . This equals the inverse 
of the original slope only in the case of perfect correlation, 
where (At/Aiy)^ = ((Aiy)2)((A/7)^). For our purposes a 
more symmetric estimate of the slope is desired, one which 
agrees with the linear regression slope in the limit of perfect 
correlation. We use simply the ratio of standard deviations 
V'((Aiy)2)/^((A/7)2) (Eq. ^). This slope was used to 
plot the dashed line in Fig. |2ja) and the full lines in Fig. [S] 
where it clearly represents the orientation of the data better.^ 
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FIG. 9: Plot of R for TIP5P water in NVT simulations with densities 
chosen to give an average pressure of one atmosphere. Not only is the 
magnitude of R low (less than 0.2) in the temperature range shown, 
but it changes sign around the density maximum. The vertical arrow 
indicates the state point used for Fig.[2](b). 



C. When competition between van der Waals and Coulomb 
interactions kills the correlation: TIP5P water 



As we shall see in the next section, the systems which show 
little correlation include several which involve both van der 
Waals and Hydrogen bonding, modeled by Lennard- Jones and 
Coulomb interactions respectively. As noted already, the lat- 
ter, being a pure inverse power-law (n = 1), by itself exhibits 
perfect correlation with slope 7 = 1/3, while the Lennard- 
Jones part has near perfect correlation. But the significant dif- 
ference in slopes means that no strong correlation is seen for 
the full interaction. To check explicitly that this is the reason 
the correlation is destroyed we have calculated the correlation 
coefficients for the Lennard- Jones and Coulomb parts sepa- 
rately in a model of water. Water is chosen because the den- 
sity of Hydrogen bonds is quite high. Simulations were done 
with the TIP5P model of watei^ which has the feature that 
the density maximum is reasonably well reproduced. This ex- 
istence of the density maximum is in fact related to pressure 
and energy becoming uncorrected, as we shall see. 

Figure [9] shows the correlation coefficients and slopes for a 
range of temperatures; the correlation is almost non-existent, 
passing through zero around where the density attains its max- 
imum value. We have separately determined the correlation 
coefficient of the Lennard- Jones part of the interaction; it 
ranges from 0.9992 at -25 °C to 0.9977 at 75 °C, even larger 
than we have seen in the SCLJ system. The reason for this 
is that the (attractive) Coulomb interaction forces the cen- 
ters of the Lennard- Jones interaction closer together than they 
would be otherwise, thus the relevant fluctuations are occur- 
ring higher up the repulsive part of the Lennard- Jones pair po- 
tential. Correspondingly the slope from this interaction ranges 
between 4.45 and 4.54, closer to the high-T, high density limit 
of 4 than was the case for the SCLJ system. This is confirmed 
by inspection of the oxygen-oxygen radial distribution func- 



tion in Ref . [26] where it can be seen that the first peak lies 
entirely to the left of the vlj = {) distance a = 0.312 nm. Fi- 
nally note that the near coincidence between the vanishing of 
the correlation coefficient and the density maximum, which 
is close to the experimental value of 4°C, is not accidental: 
The correlation coefficient is proportional to the configura- 
tional part of the thermal pressure coefficient [3v (Paper II). 
The full f3y vanishes exactly at the density maximum (4°C), 
but the absence of the kinetic term means that the correlation 
coefficient vanishes at a slightly higher temperature (~ 12° C). 



D. Results for all systems 



In Fig. 10 we summarize the results for the various sys- 
tems. Here we plot the numerator of Eq. ([6]) against the de- 
nominator, including factors of l/{kBTV) in both cases to 
make an intensive quantity with units of pressure. Since R 
cannot be greater than unity, no points can appear above the 
diagonal. Being exactly on the diagonal indicates perfect cor- 
relation {R = 1), while being significantly below it indicates 
poor correlation. Different types of symbols indicate different 
systems, as well as different densities for the same system, 
while symbols of the same type correspond to different tem- 
peratures. 

All of the simple fiquids, SCLJ, KABLJ, EXP, DB, TOL, 
show strong correlations, while METH, SPC/E, and TIP5P 
show little correlation. Values of R and 7 at selected state 
points for all systems are listed in Table [11] What determines 
the degree of correlation? The Dzugutov and TIP5P cases 
have already been discussed. The poor correlation for METH 
and SPC/E is (presumably) because these models, like TIP5P, 
involve both Lennard- Jones and Coulomb interactions. In sys- 
tems with multiple Lennard- Jones species, but without any 



O 



> 



fo 



10 



10^ 



A 

<] 

V 10 



10" 



"1 — I — I I 1 1 1 



T" 



T 1 I I I I I 



"1 1 — I I I I I I 



o 
n 
o 

A 

< 
A 

D 

O 
A 

V 
^ 

O 

♦ 
► 



SCLJ: 39.8 mol/1 

SCU: 37.4 mol/1 

SCLJ: 36.0 mol/1 

SCU: 34.6 moI/1 

SCLJ: 32.6 moW 

SCLJ Crystal: 39.8 mol/1 

EXP: 33.2 mol/1 

EXP: 30.7 mol/1 

DB 

TOL 

CU 

SQW 

KABLJ liquid 

KABLJ glass ~j 

OTP ^ 



(Forbidden region) 





A 



/ X' Increasing 
/♦ ^ Temperature 



V 



DZ 

METH 

SPC/E: 0.9 g/ml 
SPC/E: 1.0 g/ml 
SPC/E: 1.1 g/ml 
MGCU: 77 mol/1 
MGCU: 85 mol/1 
TIP5P 



* - 



* 



J \ I I I I I I 



J \ I I I I I I 



J \ I I I I I 1 1 



10" 



10" 



10 



10^ 



(<(AU)^><(AW)^>)^^^/(k TV) [GPa] 



FIG. 10: (Color online) W, U correlations for all simulated liquids; (AWAU) /{ksTV) plotted versus (((Al^)^) {{AUf)) ^'^ JOzbTV). 
Both quantities correspond to a pressure, which is given in units of GPa; for model systems not specifically corresponding to real systems, 
such as SCLJ, KABLJ, SQW, Argon units were used to set the energy and length scales. If the correlation is perfect {R — 1) the data fall on 
the diagonal marked by a dashed line. For the TIP5P model of water only temperatures with i? > are included; volumes were chosen to give 
a pressure close to zero. 



Coulomb interaction, there is overall a strong correlation be- 
cause the slope is almost independent of the parameters for a 
given kind of pair. 

As the temperature is increased, the data for the most poorly 
correlated systems, which are all hydrogen-bonding organic 
molecules, slowly approach the perfect-correlation line. This 
is consistent with the idea that this correlation is observed 
when fluctuations of both W and IJ are dominated by close 
encounters of pairs of neighboring atoms; at higher tempera- 
ture there are increasingly many such encounters, which there- 
fore come to increasingly dominate the fluctuations. Also be- 
cause the Lennard- Jones potential rises much more steeply 
than the Coulomb potential, the latter becomes less important 



as short distances dominate more. Although not obvious in 
the plot, we find that increasing the density at fixed temper- 
ature generally increases the degree of correlation, as found 
in the SCLJ case; this is also consistent with the increasing 
relevance of close encounters or collisions. 

A system quite different from the others presented so far 
is the binary square-well system, SQW, with a discontinuous 
potential combining hard-core repulsion and a narrow attrac- 
tive well (Fig. [TT} see appendix |A| for details). Such a po- 
tential models attractive colloidal systems, ^^ one of whose 
interesting features, predicted from simulations and theory, 
is the existence of two different glass phases, termed the 
"repulsive" and "attractive" glasses. ^^ The discontinuous po- 
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TABLE II: Correlation coefficients and effective slopes at selected 
state points for all investigated systems besides SCLJ. Argon units 
were used for DZ, EXP, KABLJ and SQW by choosing the length 
parameter (of the larger particle when there were two types) to be 
0.34 nm and the energy parameter to be 0.997 kJ/mol. The phase is 
indicated as liquid or glass. SQW data involves time averaging over 
periods 3.0, 3.0, 8.0 and 9.0, respectively, for the four listed state 
points. A minus sign has been included with the slope when R < 0; 
note that the 7 values only really make sense as slopes when \R\ is 
close to unity. The ensemble was NVT except for CU, DZ, MGCU, 
and SQW, where it was NVE. 



system 


p(mol/L) 


T(K) 


phase 


R 


7 


CU 


125.8 


1500 


liquid 


0.907 


4.55 


CU 


125.8 


2340 


liquid 


0.926 


4.15 


DB 


11.0 


130 


liquid 


0.964 


6.77 


DB 


9.7 


300 


liquid 


0.944 


7.45 


DZ 


37.2 


78 


liquid 


0.585 


3.61 


EXP 


30.7 


96 


liquid 


0.908 


5.98 


EXP 


33.2 


96 


liquid 


0.949 


5.56 


KABLJ 


50.7 


30 


glass 


0.858 


6.93 


KABLJ 


50.7 


70 


liquid 


0.946 


5.75 


KABLJ 


50.7 


240 


liquid 


0.995 


5.10 


METH 


31.5 


150 


liquid 


0.318 


22.53 


METH 


31.5 


600 


liquid 


0.541 


6.88 


METH 


31.5 


2000 


liquid 


0.861 


5.51 


MGCU 


85.0 


640 


liquid 


0.797 


4.74 


MGCU 


75.6 


465 


liquid 


0.622 


6.73 


OTP 


4.65 


300 


liquid 


0.913 


8.33 


OTP 


4.08 


500 


liquid 


0.884 


8.78 


OTP 


3.95 


500 


liquid 


0.910 


7.70 


SPC/E 


50.0 


200 


liquid 


0.016 


208.2 


SPC/E 


55.5 


300 


liquid 


0.065 


48.6 


SQW 


60.8 


48 


liquid 


-0.763 


-50.28 


SQW 


60.8 


79 


liquid 


-0.833 


-49.11 


SQW 


60.8 


120 


liquid 


-0.938 


-52.02 


SQW 


59.3 


120 


liquid 


-0.815 


-30.07 


TIP5P 


55.92 


273 


liquid 


-0.051 


-2.47 


TIP5P 


55.94 


285.5 


liquid 


0.000 


2.51 


TOL 


10.5 


75 
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FIG. 1 1 : Illustration of the square- well potential, indicating the four 
microscopic processes which contribute to the virial. 
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FIG. 12: (Color online) Energy-energy, virial- virial, and energy- 
virial correlation functions for SQW at packing fraction = 0.595 
and temperature T = 1.0 (normalized to unity at t = 0). The cross 
correlation has been multiplied by -1. The arrow marks the time 
t = 8, roughly 1/10 of the relaxation time (determined from the 
long-time part of the energy-virial cross-correlation function). This 
time was used for averaging. 



tential not only makes the simulations substantially different 
from a technical point of view, but there are also conceptual 
differences — in particular, the instantaneous virial is a sum of 
delta- functions in time. The dynamical algorithm employed 
in the simulations is "event-driven", where events involve a 
change in the relative velocity of a pair of particles due to 
hitting the hard-core inner wall of the potential or crossing 
the potential-step. The algorithm must detect the next event, 
advance time by the appropriate amount, and adjust the ve- 
locities of all particles appropriately. There are four kinds of 
events (illustrated in Fig. 11): (1) "collisions", involving the 



inner repulsive core; (2) "bounces", involving bouncing off 
the outer (attractive) wall of the potential; (3) "trapping", in- 



volving the separation going below the range of the outer wall 
and (4) "escapes", involving the separation increasing beyond 
the outer wall. To obtain meaningful values of the virial a 
certain amount of time- averaging must be done — we can no 
longer consider truly instantaneous quantities. As shown in 
appendix |C] the time-averaged virial may be written as the fol- 
lowing sum over all events which take place in the averaging 
interval Ly\ 
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Here r^ and Vg refer to the relative position and velocity for 
the pair of particles participating in event e, while A indicates 
the change taking place in that event. Positive contributions 
to W come from collisions; the three other event types in- 
volve the outer wall, which, as is easy to see, always gives a 
negative contribution. The default tav in the simulation was 
0.05. Strong correlations emerge only at longer averaging 
times, however. An appropriate tav may be chosen by con- 
sidering the correlation functions (auto- and cross-) for virial 



and potential energy, plotted in Fig. [T2j where the "instanta- 
neous" values E{t) and W{t) correspond to averaging over 
0.05 time units. We choose t^y c^ Tq,/ 10, where Ta is the re- 
laxation time determined from the cross-correlation function 
{AU{0)AW{t)). Data for a few state points are shown in Ta- 
ble |ll| Remarkably, this system, so different from the continu- 
ous potential systems, also exhibits strong W^ [/-correlations, 
with R = 0.94 in the case T = 1.0, ^ = 0.595 (something 



already hinted at in Fig. [12] in the fact that the curves coin- 
cide). There is a notable difference, however, compared to 
continuous systems: The correlation is negative. 

The reason for the negative correlation is that at high den- 
sity, most of the contributions to the virial are from collisions: 
a particle will collide with neighbor 1 , recoil and then collide 
with neighbor 2 before there is a chance to make a bounce 
event involving neighbor 1 . The number of collisions that oc- 
cur in a given time interval is proportional to the number of 
bound pairs which is exactly anti-correlated with the energy. 
The effective slope 7 has a large (negative) value of order -50, 
which does not seem to depend strongly on temperature. This 
example is interesting because it shows that strong pressure- 
energy correlations can appear in a wider range of systems 
that might first have been guessed. Note, however, that the 
ordinary hard-sphere system cannot display such correlations, 
since potential energy does not exist as a dynamical variable 
in this system, i.e., it is identically zero. The idea of corre- 
lations emerging when quantities are averaged over a suitable 
time interval is one we shall meet again in Paper II in the con- 
text of viscous liquids. 



IV. SUMMARY 

We have demonstrated a class of model liquids whose equi- 
librium thermal fluctuations of virial and potential energy 
are strongly correlated. We have presented detailed inves- 
tigations of the presence or absence of such correlations in 
various liquids, with extra detail presented for the standard 
single-component Lennard- Jones case. One notable aspect 
is how widespread these correlations are, appearing not just 
in Lennard- Jones potentials or potentials which closely re- 
semble the Lennard- Jones one, but also in systems involving 
many -body potentials (CU, MGCU) and discontinuous poten- 
tials (SQW). We have seen how the presence of different types 
of terms in the potential, such as Lennard- Jones and Coulomb 
interactions, spoil the correlations, even though each by itself 
would give rise to a strongly correlating system. Most simu- 
lations were carried out in the NVT ensemble; R is ensemble 
dependent, but the i? ^ 1 limit is not. 



Several of the hydrocarbon liquids studied here were sim- 
ulated using simplified "united-atom" models where groups 
such as methyl-groups or even benzene rings were represent 
by Lennard- Jones spheres. These could also be studied using 
more realistic "all-atom" models, where every atom (includ- 
ing hydrogen atoms) is included. It would be worth investi- 
gating whether the strength of the correlations is reduced by 
the associated Coulomb terms in such cases. 

In paper II we provide a detailed analysis for the single- 
component Lennard- Jones case, including consideration of 
contributions beyond the effective inverse power-law approx- 
imation and the T ^ limit of the crystal. There we also dis- 
cuss some consequences, including a direct experimental ver- 
ification of the correlations for supercritical Argon and conse- 
quences of strong pressure-energy correlations in highly vis- 
cous liquids and biomembranes. 
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APPENDIX A: DETAILS OF INTERATOMIC POTENTIALS 

Here we give more detailed information about the inter- 
atomic potentials used. These details have been published 
elsewhere as indicated, except for the case of EXP and TOL. 

CU Pure liquid Cu simulated using the many-body poten- 
tial derived from effective medium theory (EMT).E^Ell 
This is similar to the embedded atom method of Daw 
and Baskes,^^ where the energy of a given atom i, 
Ei is some nonlinear function (the "embedding func- 
tion") of the electron density due to the neighboring 
atoms. In the EMT, it is given as the energy of an 
atom in an equivalent reference system, the "effective 
medium", plus a correction term, Ei = Ec,i{ni) + 



^ref 



Y^j^i^ijiXij) -Y^jl^i^ijiXij) • Specifically, the 
reference system is chosen as an FCC crystal of the 
given kind of atom, and "equivalent" means that the 
electron density is used to choose the lattice constant 
of the crystal. The correction term is an ordinary pair 
potential involving a simple exponential, but notice that 
the corresponding sum in the reference system is sub- 
tracted (guaranteeing that the correct reference energy 
is given when the configuration is fact, the reference 
configuration). The parameters were £^0 =-3.5 10 eV; 
5o=1.413A; Vq=2A16 eV; r]2= 3.122A-^ a>:=5.178; 
A=3.602;no=0.0614A-3. 



DB Asymmetric "dumb-bell" molecules, ^^ consisting of two 
unlike Lennard- Jones spheres, labelled P and M, con- 
nected by a rigid bond. The parameters were 6^=5.726 
kJ/mol, cFp =0.4963 nm, mp=77.106 u; e^=0.66944 
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kJ/mol, (jjn =0.3910 nm, m^=15.035 u; the bond 
length was d=029 nm. Cross interactions, Cprn and 
cFpm, were set equal to the geometric and arith- 
metic means of the p and m parameters, respectively 
(Lorentz-Berthelot mixing rule23). 

DZ A monatomic liquid introduced by Dzugutov as a candi- 
date for a monatomic glass-forming system. ^^ The po- 
tential is a sum of two parts, v{r) = vi{r) + V2{r), 
with vi{r) = A{r~'^ — B) exp(c/(r — a)) for r < a 
and zero otherwise, and V2{r) = 5exp((i/(r — h)) for 
r < b, zero otherwise. The parameters are chosen to 
match the location and curvature of the Lennard- Jones 
potential: m = 16, A = 5.82, c = 1.1, a = 1.87, 
B = 1.28, 6/ = 0.27, 6 = 1.94. 

EXP A system interacting with a pair poten- 
tial with exponential repulsion U{r) = 
f [6exp(-14(r/cr - 1)) - 14(cr/r)^]. Note that 
the attractive term has the same form as the Lennard- 
Jones potential. 

KABLJ The Kob- Andersen binary Lennard- Jones liquicP^, a 
mixture of two kinds of particles A and B, with A mak- 
ing 80% of the total number. The energy and length 
parameters are caa = 1-0, cbb = 0.5, cab = 1-5, 
c^AA = 1-0, (Jbb = 0.88, aAB = 0.8. The masses 
are both equal to unity. We use the standard density 

p = i-2^ii- 

METH The Gromos 3 -site model for methanol.^ The sites 
represent the methyl (M=CH3) group (m = 15.035 
u), the O atom (m = 15.999 u), and the 0-bonded H 
atom (m = 1.008 u). The charges for Coulomb inter- 
actions are respectively 0.176 e, -0.574 e, and 0.398 e. 
The M and O groups additionally interact via Lennard- 
Jones forces, with parameters cmm = 0.9444 kJ/mol, 
eoo = 0.8496 kJ/mol, cmo = 0.9770 kJ/mol, (Jmm = 
0.3646 nm, aoo = 0.2955 nm, and (Jmo = 0.3235 
nm. Lennard- Jones interactions are smoothly cutoff be- 
tween 0.9 nm and 1.1 nm. The M-0 and 0-H distance 
is fixed at 0.136 nm and 0.1 nm, respectively, while the 
M-O-H bond angle is fixed at 108.53°. 

MGCU A model of the metallic alloy MggsCuis, simulated 
by EMT with parameters as in Ref . |22l In this poten- 
tial there are seven parameters for each element. How- 
ever, some of the Cu parameters were allowed to vary 
from their original values in the process of optimizing 
the potential for the Mg-Cu system. The parameters 
for Cu were £;o=-3.510 eV; 5o=1.413A; ^0=1.994 eV; 
r]2= 3.040A-^ K.=4.944', A=3.694; no=0.0637 A-^ 
while those for Mg were £;o=- 1-487 eV; so=1.766A; 
^0=2.230 eV; r]2= 2.54lA-i; /^=4.435; A=3.293; 
no=0.0355 A~^. It should be noted that there is an error 
in Ref. 22 : The parameter sq for Cu is given in units of 
bohr instead of A. 

OTP The Lewis-Wahnstrom three- site model of 
orthoterphenyl^^ consisting of three identical Lennard- 
Jones spheres located at the apices A B, and C of 



an isosceles triangle. Sides AB and BC are 0.4830 
nm long, while the ABC angle is 75°. The Lennard- 
Jones interaction parameters are e = 4.989 kJ/mol, 
a = 0.483 nm, while the mass of each sphere, not 
specified in Ref .^^ , was taken as one third of the mass 
of an OTP molecule, m = 76.768 u. 

SCLJ The standard single-component Lennard- Jones system 
with potential given by Eq. ( pTj ). 



SPC/E The SPC/E model of water, ^"^ in which each molecule 
consists of three rigidly bonded point masses, with an 
OH distance of 0. 1 nm and the HOH angle equal to the 
tetrahedral angle. Charges on O and each H are equal 
to -0.8476 e and +0.4238 e, respectively. O atoms inter- 
act with each other via a Lennard- Jones potential with 
6=2.601 kJ/mol and cr=0.3166 nm. 



SQW A binary model with a pair interaction consisting o f an 
infinitely hard core and an attractive square welk^^^ll] 



^ij(^) 



< dij, Vij{r) 



-uo, (Jij < r < 



(Jij{l + e), Vij{r) = 0, r > aij{l -\- e). The radius 
parameters are aAA = 1-2, (Jbb = 1, o^ab = 1-1, 
while e = 0.03 and i^o = 1. The composition was 
equimolar, and the masses of both particles were equal 
to unity. 



TIP5P In this water model^^ there are five sites associated 
with a single water molecule. One for the O atom, one 
for each H, and two to locate the centers of negative 
charge corresponding to the electron lone-pairs on the 
O. The OH bond length, and HOH bond angle are fixed 
at the gas-phase experimental values, tqh = 0.09572 
nm and 6>i:/Oi^ = 104.52°. The negative charge sites are 
located symmetrically along the lone-pair directions at 
distance tql = 0.07 nm and with an intervening angle 
Olol = 109.47°. A charge of +0.241 e is located on 
each Hydrogen site, while charges of equal magnitude 
and opposite sign are placed on the lone-pair sites. O 
atoms on different molecules interact via the Lennard- 
Jones potential with ao = 0.312 nm and eo = 0.669 
kJ/mol. 



TOL A 7-site united-atom model of toluene, consisting of 
six "ring" C atoms and a methyl group (H atoms are 
not explicitly represented). In order to handle the 
constraints more easily, only three mass points were 
used; one at the ring C attached to the methyl group 
(m = 40.065 u), and one at each of the two "meta" C 
atoms (m = 26.038) (note that with this mass distribu- 
tion, the moment of inertia is not reproduced correctly). 
Parameters were derived from the information in 
Ref.ini ering=0A602 kJ/mol, Cmethy 1=0.6694 kJ/mol, 
aring=0.315 nm, amethyi=0. 391 nm. The Lorentz- 
Berthelot rule was used for cross-interactions. ^^^ 
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APPENDIX B: CONNECTING FLUCTUATIONS TO 
THERMODYNAMIC DERIVATIVES 

If A is a dynamical quantity which depends only on the con- 
figurational degrees of freedom, then its average in the canon- 
ical ensemble (NVT) is given by (where, for convenience, we 
use a discrete-state notation, with Ai referring to the value of 
A in the ith micro-state, etc.) 






<-pU^) 



(Bl) 



where /3 = l//c^T and Q is the configurational partition func- 
tion. Then the inverse temperature derivative of {A) can be 
written in terms of equilibrium fluctuations: 



(d{A)\ Y.iAi^M-l3Ui)Ui 

\dP )y Q 

Ei Ai exp{-m) Ej exp{-pUj)U, 



(B2) 



= -({AU)-{A){U)) 
= -{AAAU). 



(B3) 
(B4) 



Now taking A = W and A = U successively we find that 



_ {AWAU) 



V dp J, 



{{Auy 



(B5) 



This last expression is precisely the formula for the slope ob- 
tained by linear-regression when plotting W against U. 

Consider now volume derivatives. Because volume depen- 
dence comes in through the micro-state values, Ai and Ui, and 
the volume derivatives of these are not necessarily related in 
a simple way, the corresponding expression is not as simple: 
The derivative of (W) with respect to volume at fixed temper- 
ature is given by 



fd{W)\ 
V dV J, 



_d_ 
W 



{{p)v- 



NUbT 



yfdM 

\dV ,^ 



(p) - Kt, 



(B6) 

(B7) 



where Kt is the isothermal bulk modulus. The derivative of 
U can be obtained by writing pressure as the derivative of 
Helmholtz free energy F (K is kinetic energy): 



ip) 



dF_ 

dV 



(^mi±m^l <Bs, 



fd{U)\ 



as 

dV 



(B9) 




FIG. 13: Illustration of replacement of discontinuous step by a finite 
slope for the square well potential for the purpose of calculating the 
virial. The limit ^ ^ is taken at the end. 



Then using 

lives of (W) and (U) 



the thermodynamic identity (^yy^ 
(3v, we obtain the ratio of volume deriva- 



KdV) 



(d{w)\ ,(m\ 

\ dV )^'\dV ), 



Kj 



(P) 



TI3v - ip) 



(BIO) 



which becomes —Kt/{TPy) when the pressure is small com- 
pared to the bulk modulus. As discussed in Paper II, Pv can 
be expressed in terms of (AUAW) again, but the fluctuation 
expression for Kt is more complicated. Thus we cannot sim- 
ply identify the lines of constant T, varying V, on a (W), {U) 
plot, as we could with lines of fixed V, varying T, by exam- 
ining the fluctuations at fixed F, T. 



APPENDIX C: VIRIAL FOR SQUARE-WELL SYSTEM 



Here we derive the expression for the time-averaged virial. 



Eq. ([17]), as a convenience for the reader. The idea is to replace 
the step uq in the potential with a a finite slope uq/S over a 
range S, and take the limit S ^ 0. We start by replacing a two- 
body interaction in three dimensions with the equivalent one- 
dimensional, one-body problem using the radial separation r 
and the reduced mass m^. Let the potential step be at r = r^ 
and define x = r — Vg (see Fig. [13]). We consider an "escape 
event" over a positive step, so that an initial (relative) velocity 
^0 becomes a final velocity vi and r goes from a value less 
than ro to a value greater than tq ^ 5. The resulting formula 
also applies for the other kinds of events. 

The contribution to the time integral of the virial from this 
event is given by 



(ro + x) 



Fdt 



ts 



(ro + x)uo 
3S 



dt 



(CI) 



where F is the (constant) force in the region < x < S and 
ts is the time taken for the 'particle' (the radial separation) to 
cross this region. The trajectory x{t) is given by the standard 
formula for uniform acceleration 
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1 1^0 2 

X{t) =vot- 1 , 

2 drrir 



(C2) 



which by setting x{ts) = S gives the following expression for 

t6- 



-rrirV, 



1 



r^^O 



-rrirVi + uq 



from which Av is obtained as 



(C6) 



t5 



rrirVo 

Uo 



rrirVo 



2mr 

Uo 



(C3) 



here we have taken the negative root, appropriate for a positive 
Uq (we want the smallest positive ts). Returning to A, it can 
be split into two parts as follows 



Av = vi - ^0 = \/vo^ - 2uo/mr - vo, (C7) 



thus the expression for A becomes 



roTTlr 



Av 



rrir 



r • Av 



(C8) 



-fr^--f: 



x{t)uo 
35 



dt 



(C4) 



Consider the second term: using the expression for x{t) 
from Eq. (|C2|), we see that the result of the integral will in- 



volve a term proportional to t^ and one proportional to t|. 



Using Eq. ( |C3[ ) to replace ts ex 5, and noting the 5 in the 
denominator, the terms will have linear and quadratic depen- 
dence on 6, respectively. Thus they will vanish in the limit 
6^0. On the other hand, the first term gives 



rpup 
' 36 



ts 




(C5) 



This expression can be simplified by writing it in terms of 
the change of velocity Av = Vi — V{).lx\ the one-body problem 
conservation of momentum does not hold, and vi is given by 
energy conservation: 



where in the last expression a switch to three-dimensional no- 
tation was made, recognizing that for central potentials Av 
will be parallel to the displacement vector between the two 
particles. This expression, derived for escape events, must 
also hold for capture events since these are time-reverses of 
each other, and the virial is fundamentally a configurational 
quantity, independent of dynamics (the above expression is 
time-reversal invariant because the change in the radial com- 
ponent of velocity is the same either way, since although the 
"initial" and "final" velocities are swapped, they also have op- 
posite sign). Bounce and collision events may be treated by 
dividing the event into two parts at the turning point (where 
the relative velocity is zero), noting that each may be treated 
exactly as above, then adding the results back together. If we 
now consider all events that take place during an averaging 
time tav), we get the time-averaged virial as 



W 



otp. 



E 



m 



r,e^e ' ■^'^e 



(C9) 
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